VMI class demo

23/11/20 PH

The VMI class builds on the tmoDataBase class for VMI data analysis.

Setup

In [1]:
import numpy as np

from h5py import File

from pathlib import Path

# HV imports
import holoviews as hv
from holoviews import opts
hv.extension('bokeh', 'matplotlib')

import xarray as xr


# Memory profiler - OPTIONAL for testing
# https://timothymonteath.com/articles/monitoring_memory_usage/
%load_ext memory_profiler
%memit


# Quick hack for slow internet connection!
%autosave 36000


# Import class - this requires pointing to the `tmoDataBase.py` file.
# See https://github.com/phockett/tmo-dev

# Import class from file
import sys
sys.path.append(r'/cds/home/p/phockett/dev/tmo-dev/tmo')
import tmoDataBase as tb
from vmi import VMI  # VMI class

tb.setPlotDefaults(fSize=[500,500])  # Set plot defaults (optional)
peak memory: 162.18 MiB, increment: 0.30 MiB
Autosaving every 36000 seconds
In [2]:
# Set main data path
baseDir = Path('/reg/data/ana15/tmo/tmolw0618/scratch/preproc/v2')

# Setup class object and which runs to read.
data = VMI(fileBase = baseDir, runList = range(89,97+1), fileSchema='_preproc_elecv2')

# The class has a bunch of dictionaries holding info on the data
# data.runs['files']

# Read in the files with h5py
# There is very basic checking implemented currently, just to see if 'energies' is present.
data.readFiles()
*** WARNING: key 90 missing energies data, will be skipped.
*** WARNING: key 92 missing energies data, will be skipped.
*** WARNING: key 94 missing energies data, will be skipped.
*** WARNING: key 95 missing energies data, will be skipped.
*** WARNING: key 96 missing energies data, will be skipped.
Read 9 files.
Good datasets: [89, 91, 93, 97]
Invalid datasets: [90, 92, 94, 95, 96]

Generate VMI data

  • This uses np.histogram2d^ to construct images from two passed data dimensions, the default is the electron data ['yc','xc'].
  • Multiple filters can be set,

^ Faster methods to be developed, this is currently a bit slow (10s of seconds).

Filters

  • The defaults as gas on and gas off as signal and background filtersets.
  • These are currently set as a nested dictionary in self.filters, and can be set manually, or with self.setFilters() (additive, or reset all).
  • These are currently set as matched values, or inclusive ranges.

More filter options to come!

In [3]:
# Default filter settings
data.filters
Out[3]:
{'signal': {'gas': True, 'desc': 'Signal filter.'},
 'bg': {'gas': False, 'desc': 'Background filter.'}}
In [4]:
# Add a filter to an existing filterset
data.setFilter({'signal':{'energies':[0.05, 0.1]}})
data.filters
Out[4]:
{'signal': {'gas': True, 'desc': 'Signal filter.', 'energies': [0.05, 0.1]},
 'bg': {'gas': False, 'desc': 'Background filter.'}}
In [5]:
# Add a filter to all filtersets
data.setFilter({'ts':3.185})
data.filters
Out[5]:
{'signal': {'gas': True,
  'desc': 'Signal filter.',
  'energies': [0.05, 0.1],
  'ts': 3.185},
 'bg': {'gas': False, 'desc': 'Background filter.', 'ts': 3.185}}
In [6]:
# Add a new filterset
data.setFilter({'ref':{'gas':True, 'energies':0.05,'ts':3.185}})
data.filters
Out[6]:
{'signal': {'gas': True,
  'desc': 'Signal filter.',
  'energies': [0.05, 0.1],
  'ts': 3.185},
 'bg': {'gas': False, 'desc': 'Background filter.', 'ts': 3.185},
 'ref': {'gas': True, 'energies': 0.05, 'ts': 3.185}}
In [7]:
# For a full filter reset pass reset=True
# This is currently set to just ckear all filters
data.setFilter(reset=True)
data.filters
Out[7]:
{}
In [8]:
# Manual setting
data.filters = {'signal':{'gas':True,
                          'desc': 'Signal filter.'},
                'bg':{'gas':False,
                      'desc': 'Background filter.'},
                }

Run image generation

This defaults to all data and all filtersets, so may take a little while (~seconds to minutes).

In [9]:
# dims = ['yc','xc']  # Optionally pass dimensions to use from the raw data, the default is ['yc','xc']
data.genVMIXmulti()
Generating VMI images for filters: signal
Generating VMI images for filters: bg

This will generate images and output to an Xarray Dataset, self.imgStack, with dimensions (yc,xc,run) for each filterset.

In [10]:
data.imgStack
Out[10]:
<xarray.Dataset>
Dimensions:  (run: 4, xc: 1048, yc: 1048)
Coordinates:
  * yc       (yc) float64 -0.5 0.5 1.5 2.5 ... 1.044e+03 1.046e+03 1.046e+03
  * xc       (xc) float64 1.046e+03 1.046e+03 1.044e+03 ... 1.5 0.5 -0.5
  * run      (run) int64 89 91 93 97
    norm     (run) int64 72318 73364 71472 36198
Data variables:
    signal   (yc, xc, run) float64 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    bg       (yc, xc, run) float64 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0

Filter parameters and image generation settings and metrics are set to self.data[run][filterset].

In [11]:
data.data[89]['signal']
Out[11]:
{'metrics': {'shots': (72318, 1000),
  'selected': 61985,
  'gas': 61985,
  'events': 61985000,
  'norm': 72318},
 'mask': array([ True, False,  True, ..., False, False,  True])}

Image manipulation

Computation

Operations can be carried out directly on the imgStack. Note that, by default, the images are normalised to number of shots, but not by any other parameters.

In [12]:
# Subtraction, and set as new image set
data.imgStack['sub'] = data.imgStack['signal'] - data.imgStack['bg']
data.imgStack
Out[12]:
<xarray.Dataset>
Dimensions:  (run: 4, xc: 1048, yc: 1048)
Coordinates:
  * yc       (yc) float64 -0.5 0.5 1.5 2.5 ... 1.044e+03 1.046e+03 1.046e+03
  * xc       (xc) float64 1.046e+03 1.046e+03 1.044e+03 ... 1.5 0.5 -0.5
  * run      (run) int64 89 91 93 97
    norm     (run) int64 72318 73364 71472 36198
Data variables:
    signal   (yc, xc, run) float64 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    bg       (yc, xc, run) float64 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    sub      (yc, xc, run) float64 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0

(More to follow, but see the Xarray docs for general info.)

In [13]:
# The data can also be output as a 4D array, (x,y,run,name), which may be more useful for further computations
data.restackVMIdataset()
Out[13]:
<xarray.DataArray 'stacked' (yc: 1048, xc: 1048, run: 4, name: 3)>
array([[[[0., 0., 0.],
         [0., 0., 0.],
         [0., 0., 0.],
         [0., 0., 0.]],

        [[0., 0., 0.],
         [0., 0., 0.],
         [0., 0., 0.],
         [0., 0., 0.]],

        [[0., 0., 0.],
         [0., 0., 0.],
         [0., 0., 0.],
         [0., 0., 0.]],

        ...,

        [[0., 0., 0.],
         [0., 0., 0.],
         [0., 0., 0.],
...
         [0., 0., 0.],
         [0., 0., 0.],
         [0., 0., 0.]],

        ...,

        [[0., 0., 0.],
         [0., 0., 0.],
         [0., 0., 0.],
         [0., 0., 0.]],

        [[0., 0., 0.],
         [0., 0., 0.],
         [0., 0., 0.],
         [0., 0., 0.]],

        [[0., 0., 0.],
         [0., 0., 0.],
         [0., 0., 0.],
         [0., 0., 0.]]]])
Coordinates:
  * yc       (yc) float64 -0.5 0.5 1.5 2.5 ... 1.044e+03 1.046e+03 1.046e+03
  * xc       (xc) float64 1.046e+03 1.046e+03 1.044e+03 ... 1.5 0.5 -0.5
  * run      (run) int64 89 91 93 97
    norm     (run) int64 72318 73364 71472 36198
  * name     (name) <U6 'signal' 'bg' 'sub'

Visualisation

Holoviews plus Bokeh or Plotly backends can be used to quickly construct interactive visualisation (including interactive HTML versions, although file sizes may become large).

In [14]:
# To reduce network bandwidth requirements, for those with slow internet connections, the image stack can be downsampled prior to visualisation
data.downsample(step=[5,5])

# This is set to self.imgReduce, note this preserves coords.
data.imgReduce
Out[14]:
<xarray.Dataset>
Dimensions:  (run: 4, xc: 209, yc: 209)
Coordinates:
  * yc       (yc) float64 1.5 6.5 11.5 16.5 ... 1.032e+03 1.036e+03 1.042e+03
  * xc       (xc) float64 1.044e+03 1.04e+03 1.034e+03 1.03e+03 ... 14.5 9.5 4.5
  * run      (run) int64 89 91 93 97
    norm     (run) int64 72318 73364 71472 36198
Data variables:
    signal   (yc, xc, run) float64 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    bg       (yc, xc, run) float64 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    sub      (yc, xc, run) float64 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
In [15]:
# Single selected image by run and type
# Use the interactive histogram to adjust the colormapping, or pass clims = [min, max] to set explicitly
# Use the menu bar to pan, zoom, save and reset the image.
data.showImg(run=89, name='signal')
WARNING:param.RasterPlot04576: Due to internal constraints, when aspect and width/height is set, the bokeh backend uses those values as frame_width/frame_height instead. This ensures the aspect is respected, but means that the plot might be slightly larger than anticipated. Set the frame_width/frame_height explicitly to suppress this warning.
In [16]:
# For further use, the hv.Image object created can be returned and used...
# (See the HoloViews docs for more info: https://holoviews.org/getting_started/Introduction.html)
imgPlot = data.showImg(run=89, name='signal', returnImg=True)

# Plot image + 1D slices
# Note here the slices are selected by coods, and the plot axes are linked.
(imgPlot.hist() + imgPlot.sample(xc=500) + imgPlot.sample(yc=552)).cols(1)
WARNING:param.RasterPlot04977: Due to internal constraints, when aspect and width/height is set, the bokeh backend uses those values as frame_width/frame_height instead. This ensures the aspect is respected, but means that the plot might be slightly larger than anticipated. Set the frame_width/frame_height explicitly to suppress this warning.
Out[16]:
In [17]:
# All images can be viewed with a HoloMap + widgets...

# # Firstly set to an hv.Dataset
# imgDS = hv.Dataset(data.restackVMIdataset())

# # Then a HoloMap of images
# imgDS.to(hv.Image, kdims=['xc','yc']).opts(colorbar=True).hist()

data.showImgSet()
WARNING:param.RasterPlot09716: Due to internal constraints, when aspect and width/height is set, the bokeh backend uses those values as frame_width/frame_height instead. This ensures the aspect is respected, but means that the plot might be slightly larger than anticipated. Set the frame_width/frame_height explicitly to suppress this warning.
dir(data)